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Abstract. We study the phase diagram of the frustrated t—t' Hubbard model on the 
square lattice by using a novel variational wave function. Taking the clue from the backflow 
correlations that have been introduced long-time ago by Feynman and Cohen and have been 
used for describing various interacting systems on the continuum (like liquid '^He, the electron 
jellium, and metallic Hydrogen), we consider many-body correlations to construct a suitable 
approximation for the ground state of this correlated model on the lattice. In this way, a very 
accurate ansatz can be achieved both at weak and strong coupling. We present the evidence 
that an insulating and non-magnetic phase can be stabilized at strong coupling and sufficiently 
large frustrating ratio t' /t. 



The Hubbard model on the square lattice with nearest- and next-nearest-neighbor hoppings 
has been widely studied by many authors with different numerical techniques and contradictory 
outcomes. The Hamiltonian is given by: 



where q^(ci^o-) creates (destroys) an electron with spin a on site i, rii^a- = cJ^q^o-, Uj is 
the hopping amplitude (denoted by t and t' for nearest- and next-nearest-neighbor sites, 
respectively), and U is the on-site Coulomb repulsion. In the following, we will consider the 
half-filled case, where the number of electrons is equal to the number of sites. The model 
of Eq. ([1]) represents a simple prototype for frustrated itinerant materials. In the presence of a 
finite t' /t, there is no more a perfect nesting condition that leads to antiferromagnetism for any 
finite U and non-conventional phases may be stabilized at zero temperature, like for instance 
spin liquids with no magnetic order. 

The first numerical study of this model is due to Lin and Hirsch, [Ij who found the existence 
of a critical Uc for the appearance of antiferromagnetism at finite values of t' /t. More recent 
studies have been done by Imada et al., [21 |3l Hj by using the Path Integral Renormalization 
Group approach, by Yokoyama et ai, [5] by using a variational Monte Carlo method, and by 
Tremblay et al., [6] by a Variational Cluster Approximation. Remarkably, all these numerical 
approaches give very different results for the ground-state properties of this simple correlated 
model. In fact, there are huge discrepancies for determining the boundaries of various phases, 
but also for characterizing the most interesting non-magnetic insulator. Furthermore, also the 
possibility to have superconductivity at small values of U/t is controversial. 





In the following, we will show our numerical results, which are based upon an improved 
variational Monte Carlo approach that contains backflow correlations. [7] Before doing that, 
it is useful to remind how to construct suitable variational states to describe different phases. 
Variational wave functions for the unfrustrated Hubbard model, which has antiferromagnetic 
long-range order, can be easily constructed by considering the ground state \AF) of a mean-field 
Hamiltonian containing a band contribution and a magnetic term 

nAF = -Yl + B.C. + AaF E "^^^ (2) 

where is the x component of the spin operator Sj = {Sj,Sj,Sj) and Q is a suitable 
pitch vector, e.g., Q = (vr,7r) for the Neel phase. In order to have the correct spin-spin 
correlations at large distance, we have to apply a suitable long-range spin Jastrow factor, namely 
\^af) = J^s\AF), with Js = ex.p[—^J2i.j'^i,jSiSj]j which governs spin fluctuations orthogonal 
to the magnetic field A^p. [8j It is important to stress that the mean-field state \AF) can easily 
satisfy the single-occupancy constraint by taking A^f — > oo; in this limit, it also contains the 
virtual hopping processes, which are generated by the kinetic term, implying that it is possible 
to reproduce super-exchange processes. 

On the other hand, in pure spin models, namely when U is infinite and charge fluctuations 
are completely frozen, spin-liquid (i.e., non- magnetic) states can be constructed by considering 
the ground state \BCS) of a BCS Hamiltonian with singlet pairing 

hbcs = - e tijia^^j,- + H.c. + E A'icsiciAi + 4A,0 + (3) 

and then applying to it the so-called Gutzwiller projector, \RVB) = Vg\BCS), where Vg = 
Yliil — grii^^ni^i) and 5 = 1. [9] These kind of states can be remarkably accurate and represent 
important tools for the characterization of disordered spin-liquid ground states. [lOllllj However, 
whenever U jt is finite, the variational state must also contain charge fluctuations. In this 
regard, the simplest generalization of the Gutzwiller projector with 5 7^ 1, which allows doubly 
occupied sites, is known to lead to a metallic phase. [12] One particularly simple way to obtain 
a Mott insulator with no magnetic order is to add a sufficiently long-range Jastrow factor 
J = exp[— i X^ij ■^ij^i'^i]) = So- ^«,o- being the local density. [13] Nevertheless, the accuracy 
of the resulting wave function \^bcs) = J\BCS) can be rather poor in two dimensions for large 
on-site interactions, especially in presence of frustration, since the super-exchange energy scale 
is not correctly reproduced. In fact, in contrast to the previous case with magnetic order, within 
the uncorrelated state \BCS) it is not possible to avoid a finite amount of double occupancies, 
and the Gutzwiller factor is mandatory to project out high-energy configurations. Here, we 
propose a simple improvement of (general) correlated wave functions in order to mimic the 
effect of virtual hoppings, leading to the super-exchange mechanism. In particular, we want 
to modify the single-particle orbitals, in the same spirit of backflow correlations, which have 
been proposed long time ago by Feynman and Cohen to obtain a quantitative description of the 
roton excitation in liquid Helium. yL4j The backflow term has been implemented within quantum 
Monte Carlo calculations to study bulk liquid ^He, [151 IS] a-nd used to improve the description 
of the electron jellium both in two and three dimensions. |17^ ll8j More recently, it has been also 
applied to metallic Hydrogen. [19j 

Originally, the backflow term corresponds to consider flctitious coordinates of the electrons 
r^, which depend upon the positions of the other particles, so to create a return flow of current: 



= ra + E ^a,/3N (r/3 " ^a) , 

/3 



(4) 



where are the actual electronic positions and r/Q,/3[x] are variational parameters depending 
in principle on all the electronic coordinates, namely on the many-body configuration The 
variational wave function is then constructed by means of the orbitals calculated in the new 
positions, i.e., 0(r^). Alternatively, the backfiow term can be introduced by considering a linear 
expansion of each single-particle orbital: 

M^i) - Ki^a) = 0fc(ra) + J2 Cq,/3N M^Is), (5) 

where Cq,^/3[x] are suitable coefficients that depend on all electron coordinates. The definition ([5]) 
is particularly useful in lattice models, where the coordinates of particles may assume only 
discrete values. In the Hubbard model, the form of the new "orbitals" can be fixed by considering 
the U ^ t limit, so to favor a recombination of neighboring charge fluctuations (i.e., empty and 
doubly-occupied sites): 

(f)k{ri,a) = ^4>k{ri,a) + r]^tijDiHj(l)k{rj^„), (6) 

j 

where we used the notation that 4>k{^i^a) = (0|ci^(j|</>A:), \4>k) being the eigenstates of the mean- 
field Hamiltonian ([2]) or ([3]), Di = Ui^i^rii^i, Hi = hi^^hi^i, with hi^a- = (1 — n'i,a), so that Di and 
Hi are non zero only if the site i is doubly occupied or empty, respectively; finally e and rj are 
variational parameters (we can assume that e = 1 if DiHj = 0). As a consequence, already 
the determinant part of the wave function includes correlation effects, due to the presence of 
the many body operator DiHj. The previous definition of the backfiow term preserves the spin 
SU(2) symmetry. A further generalization of the new "orbitals" can be made, by taking all the 
possible virtual hoppings of the electrons: 

(piiri^^) = e0fc(ri,^) + r]i^tijDiHj(j)kirj,a) + 

j 

m X! tijni,ahi-anj-ahj^^<pk{rj^a) + % X! % {Dirij-crhj^a + Ui^ahi-aHj) 4>k{rj^^), (7) 

j j 

where e, rji, r]2, and r/3 are variational parameters. The latter two variational parameters are 
particularly important for the metallic phase at small U/t, whereas they give only a slight 
improvement of the variational wave function in the insulator at strong coupling. For simplicity, 
we take the same parameters for up and down electrons. 

Thanks to backfiow correlations, it is possible to obtain a correct extrapolation to the infinite- 
U limit (i.e., to the variational energy obtained with the fully projected states g' = 1 in the 
Heisenberg model). On the contrary, without using backfiow terms, the energy of the BCS 
state, even in presence of a fully optimized Jastrow factor, is few hundredths of J = 4t^/C7 
higher than the expected value, see Fig. [H The importance of backfiow correlations is even 
more evident in the frustrated case, where they are essential also for improving the accuracy of 
the antiferromagnetic wave function. 

In order to draw the ground-state phase diagram of the t—t' Hubbard model, we consider 
three different wave functions, all with backfiow correlations: One non-magnetic state \^bcs) 
and two antiferromagnetic states I^'af) with pitch vectors Q = (vr, vr) and Q = (vr,0), relevant 
for small and large t' /t, respectively. The variational phase diagram is reported in Fig. [5J 

The important outcome is that without backfiow terms, the energies of the spin-liquid wave 
function are always higher than those of the magnetically ordered states, for any value of the 
frustration t' /t. Instead, by inserting backfiow correlations, a spin-liquid phase can be stabilized 
at large enough U jt and frustration. For example, this can be seen in Fig. [3l where we show 




Figure 1. Variational energies per site (in unit of J = At^ /U) for the BCS state with the 
Jastrow factor, with and without backflow correlations, for 98 sites and t' = (left panel) and 
t' /t = 0.7 (right panel). The results for the wave function with antiferromagnetic order and no 
BCS pairing are also shown. Arrows indicate the variational results obtained by applying the 
full Gutzwiller projector (i.e., 5 = 1) to the mean-field states for the corresponding Heisenberg 
models. 
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Figure 2. Left panel: Phase diagram as obtained by comparing the variational energies of 
different wave functions, all with backflow correlations. Right panel: Phase diagram as obtained 
by Yokoyama et al. by their variational approach. |5j 



the variational energies for the three aforementioned wave functions with and without backflow 
correlations, at U /t = 16. 

In order to study the metal-insulator transition, we look at the static density-density 
correlations N{q) = {n-gUq) (where Ug is the Fourier transform of the local density n^). Indeed, 
N{q) shows a linear behavior for \q\ ^ in the metallic phase and a quadratic behavior in 
the insulating region. [T3] For small Coulomb repulsion and finite t'/t, N{q) has the linear 
behavior for |g| — > 0, typical of a conducting phase. Further, a very small superconducting 
parameter with dx2-y2 symmetry can be stabilized (e.g., A^^^ = ztA^c'^ at nearest neighbors) 
suggesting that long-range pairing correlations, if any, are tiny. In this respect, we compare 
in Table [T] the optimized Abcs when the spin- liquid wave function \^bcs) contains or not 
backflow correlations, for various U /t and t'/t = 0.75. Data show that when accuracy increases, 
by means of backflow correlations, the BCS pairing is reduced by an order of magnitude. By 
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Figure 3. Variational energies for three different wave functions on the frustrated square lattice 
with and without backflow correlations (left and right panels, respectively). Red dots denote 
the energies of the spin-liquid wave function, blue triangles the energies of the magnetic state 
with Q = (7r,7r), and green squares the energies of the magnetic state with Q = (vr, 0). Data 
are shown for U/t = 16 and 98 sites. 
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Figure 4. Variational results for N{q) divided by \q\ for 98 (empty symbols) and 162 (full 
symbols) sites and t'/t = 0.75. The metal-insulator transition takes place between U/t = 7 and 
U/t = 8, where N{q) changes from a linear to quadratic behavior for \q\ 



increasing U/t, a metal-insulator transition is found and N{q) acquires a quadratic behavior in 
the insulating phase, indicating a vanishing compressibility. In Fig. HI we show the variational 
results for N{q) as a function of U/t for t'/t = 0.75. The insulator just above the transition is 
magnetically ordered and the variational wave function has a large Aaf', the transition is likely 
to be first order, since the parameter A^p has a jump across the metal- insulator transition. 

In the frustrated regime with t'/t ~ 0.7, by further increasing U/t, there is a second transition 
to a disordered insulator. Indeed, for U/t > 14, the energy of the BCS wave function becomes 
lower than the one of the antiferromagnetic state. In this respect, the key ingredient to have 
such an insulating behavior is the presence of a long-range Jastrow term J', which turns a 
BCS superconductor into a Mott insulator. [TS] It should be noted that the spin liquid wave 
function contains a superconducting gap with d^2_y2 symmetry, in contrast to what was found 
in the infinite-f7 limit, namely in the frustrated Heisenberg model, by a similar variational 
approach. [lOj In fact, in the latter case, the BCS parameter contains both a term with d^2_y2 
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Table 1. BCS pairing Abcs for various U/t in the metalhc region at t' /t = 0.75 and 98 sites. 



Wave function 


Energy {U/t = 2Q,t' /t = 0.7) 


Energy {U/t = 8, t' /t = 0.3) 


J\BCS) 


-0.1950(1) 


-0.4016(1) 


JhdJ\BCS) 


-0.2061(1) 


-0.4180(1) 


J\BCS + Backflow) 


-0.23516(4) 


-0.4879(1) 


+ Backflow) 


-0.23257(3) 


-0.5222(1) 



Table 2. Variational energies (in unit of t) for three spin-liquid wave functions and for the best 
antiferromagnetic state with Neel order. The cluster contains 98 sites. Jhd is a short-range 
many-body Jastrow factor that has been used in Ref. [5]. 



symmetry and a further term with dxy symmetry. However, the energy gain due to the latter 
term is very small in the Heisenberg model (i.e., order of O.OOIJ) and it is very hard to detect it 
when charge fluctuations are allowed (i.e., in the Hubbard model). In all cases that have been 
analysed, we found that the dxy term is not stable in the thermodynamic limit, but converges 
to zero as the number of lattice sites is increased. 

We can make a direct comparison of our energies with the ones obtained by Yokoyama and 
collaborators, [5j who used a similar variational wave function containing a particular many-body 
Jastrow factor Jhd to correlate empty and doubly occupied sites at nearest-neighbor distances. 
In TableO we report the variational energy of the simple spin-liquid state \^ bcs) = J\BCS), 
together with the improved energies, which are obtained by adding the Jastrow term Jhd or 
by considering backflow correlations. We notice that the latter state always gives much lower 
energies than the one obtained with the additional Jastrow factor. In particular, let us consider 
the case of U/t = 8 and t'/t = 0.3, which should be magnetically ordered according to our 
calculations and disordered according to Ref. [5] (see Fig. [2|). In this case, even though our spin- 
liquid state has a much better energy than the one with Jhd, the best wave function has Neel 
order, indicating that the stability region of antiferromagnetism is larger than what predicted 
by Yokoyama and collaborators. [5] 

To conclude, we compare the variational energies with the ones obtained by the Green's 
function Monte Carlo approach, implemented within the Fixed Node (FN) approximation. In 
brief, the FN method allows one to filter out the high-energy components of a given state and 
to find the best variational state with the same nodes of the starting one. [20] On the lattice, 
the FN method can be defined as follows: Starting from the original Hamiltonian 7Y, we define 
an effective Hamiltonian by adding a perturbation O: 

n^ff = n + o. (8) 

The operator O is defined through its matrix elements and depends upon a given guiding function 
1^'), that is for instance the variational state itself 



_ j —T~(-x',x if Sx'^x = '^x'T~(-x',x'^x > 

[ J2y,sy,x>o'^y,x-§^ for X = X, 



(9) 
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Figure 5. Comparison between the variational (VMC) energies per site (with backflow 
correlations) and the FN ones. Data are shown for U /t = 16 and 98 sites. 



where = (^^l^)) 1^;) being a generic many-body configuration. The most important property 
of this effective Hamiltonian is that its ground state |^'o) can be efficiently computed by 
using the Green's function Monte Carlo technique, |21l which allows one to sample the 
distribution n^,. oc {x\^){x\^q) by means of a statistical implementation of the power method: 
n oc lim„_^oo G"!!'^, where 11° is a starting distribution and G^'^x = '^x'{J^^x',x — is 
the so-called Green's function, defined with a large or even infinite positive constant A, 5xi ^x 
being the Kronecker symbol. The statistical method is very efficient since in this case all the 
matrix elements of G are non-negative and, therefore, it can represent a transition probability 
in configuration space, apart for a normalization factor hx = X^x' Gx\x- In this case, it follows 
immediately that the asymptotic distribution 11 is also positive and, therefore, we arrive at the 
important conclusion that the ground state of Ti^^^ has the same signs of the chosen guiding 
function (i.e., the best variational state). 

In Fig. m we show the variational energies per site (with backfiow correlations) and the FN 
ones for C//t = 16 on a 98-site lattice. The small energy difference between the pure variational 
energies and the FN ones demonstrates the accuracy of the backfiow states. Notice that I^'af) 
and \^ Bcs) have different nodal surfaces, implying different FN energies. 

In order to verify the magnetic properties obtained within the variational approach, we can 
consider the static spin-spin correlations S{q) = (SgS^g), where is the Fourier transform of 
the local spin Sf. Although the FN approach may break the SU(2) spin symmetry, favoring a 
spin alignment along the z axis, S{q) is particularly simple to evaluate within this approach, [20] 
and it gives important insights into the magnetic properties of the ground state. In Fig. EJ we 
report the comparison between the variational and the FN results by using the non-magnetic 
state \^BCs)- Remarkably, in the unfrustrated case, where antiferromagnetic order takes place, 
the FN approach is able to increase spin-spin correlations at Q = {tt,tt), even by considering 
the non-magnetic wave function to fix the nodes. In this case, the FN results are qualitatively 
different from the pure variational ones, which indicate no magnetic order in the thermodynamic 
limit. A finite value of the magnetization is also plausible in the insulating region just above the 
metallic phase at strong frustration (i.e., t' /t ~ 0.75), confirming our variational calculations. On 
the contrary, by increasing electron correlation, FN results change only slightly the variational 
value of S'(7r,7r), indicating the stability of the disordered state. Therefore, the FN results 
confirm that a spin liquid region can be stabilized only at large enough U/t, while the insulator 
close to the metallic region is magnetically ordered. 

In summary, the backflow wave functions represent simple and useful generalizations of 
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Figure 6. Variational (empty symbols) and FN (full symbols) results for ^(vr, tt)/N, for N = 18, 
50, 98, and 162. All the calculations have been done by using the projected BCS wave function; 
U/t = 16 and t'/t = (triangles), U/t = 24 and t'/t = 0.7 (circles) and U/t = 8 and t'/t = 0.75 
(squares). Lines are guides to the eye. 



standard projected states, which are used to describe strongly correlated materials. They are 
highly accurate and may give important insights into the ground-state properties of frustrated 
models with itinerant electrons. 
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